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ABSTRACT 

The far-field radiation pattern of an antenna may be determined as the fast Fourier transform 
(FFT) of the aperture distribution. When the antenna is electrically large and a detailed pattern in 
(wo dimensions is required, computer run-times exceeding an hour may result. A more tinie- 
efficient algorithm for computing the discrete Fourier transform (DFT) may result in a more effec- 
tive analysis and design process. Significant savings in epu time will improve the computer turn- 
around time and circumvent the need to resort to weekend runs. 

A FORTRAN program to calculate the DFT using the Winograd Fourier transform algorithm 
was adapted to the IBM 360/91 computer and extended to handle complex input data vectors up to 
length N = 5040, Transforms may be computed for any data length given by the product of four 
mutually prime numbers selected from the integers 16, 9, 8, 7, 5, 4, 3, and 2 (e.g., N = 9 , 8*7*5 = 
2520). The WFT was used to compute antenna patterns for cophase and linear phase gradient 
apertures. The results were essentially identical to those previously computed with a conventional 
radix-2 FFT. 

Significant lime savings were realized with the WFT program. Run-time comparisons were 
made between WFT lengths of 1008, 2520 and 5040 and FFT lengths of 1024, 2048, and 4096, 
respectively. A minimum 4.6 to 1 speed advantage was demonstrated over this range. On the basis 
of the WFT timings it was estimated that two-dimensional transforms would require about one min- 
ute, ten minutes, and 40 minutes for the 1008, 2520, and 5040 point transforms, respectively. This 
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is sufficient to make two-dimensional transforms up to N = 2520 feasible within reasonable com- 
puter run-time limitations. 
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FAR-FIELD RADIATION PATTERNS OF APERTURE ANTENNAS 
BY THE WINOGRAD FOURIER TRANSFORM ALGORITHM 


INTRODUCTION 

The computation offar-field radiation patterns of aperture antennas using the Fourier trans- 
form is well documented in the literature (Ref, 1 and 2) and in current use. The fast Fourier trans- 
form (FFT) algorithm is used in this application because of its computational speed advantage over 
the conventional discrete Fourier transform (DFT). Even so, computer run times exceeding one 
hour are expected as analysis techniques arc extended to two-dimensional problems. A significant 
reduction in epu run time would allow more effective analysis and design as computer runs may not 
be so time consuming as to be relegated to weekend runs only. Furthermore, the financial savings 
may be significant. 

Recent studies (Ref. 3 and 4) have explored the feasibility of using the fast Walsh transform 
(FWT) for this application. While time savings on the order of ten to one were reported, the FWT 
gave only marginally satisfactory results for real data and failed to produce the normal squint associ- 
ated with a linear phase gradient on the aperture distribution. Additionally, there were serious diffi- 
culties in calibrating the scqucncy axis, 

Other investigators (Ref. 5, 6, 7, 8, 9, and 10) have reported on a new algorithm for computing 
the DFT. The Winograd Fourier transform algorithm (WFT) requires substantially fewer multiplica- 
tions than the FFT while the number of additions, for some cases, remains near FFT levels. The 
WFT may be used effectively on short data sequences where the number of elements is prime. For 
long transform lengths, however, the number of additions becomes excessive and a direct applica- 
tion of the WFT becomes impractical. In this case it is useful to employ a combination of “smali-N”/ , 
WFT algorithms with a multidimensional expansion to extend the range of application to data 
lengths in excess of several thousand. 

FROM DISCRETE FOURIER TRANSFORM TO DISCRETE CONVOLUTION 

The DFT of a data vector x(n) = jx(0), x(l), • • • • x(N- 1)| is defined as: 
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k = 0, 1 , 2, ■ • • N - l 
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Tlic DPT is to discrete-time signals what the Fourier transform 
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is to continuous-time signals. The extension of equation (2) to the DFT of equation (1) ci'n be 
viewed intuitively. 

liquation (I) describes the generation of N equations from which the elements in the DFT 
X(0), X(l), • • * X(N- I) may be computed: 
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Using the FFT algorithm to evaluate this matrix product results in a very substantial computational 
reduction over a direct calculation. The complexity of computation arises from the (N - 1) by 
(N - 1) lower right-hand section of the transform matrix: 
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N-! 

X(k) = x(n)w kn k = 1 , 2, • • • N - 1 
n=l 

From X(k) we can retrieve X(k) by 

N-l 

X(0) = x(n) 

li n 0 

X(k) « x(0) + 2(k) k= t, 2, • • ' N-l 


(4) 


(5) 


Consider the case where N = 5. Equation (4) then becomes: 
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where the exponents of w are written modulo 5, i,c., w® = w^ ■ w 3 = w 3 , as w** = 1 . A simple per- 
mutation of rows and columns in (6) will demonstrate how the DFT process may be converted to 


cyclic convolution. First interchange the last two columns and then the last two rows: 
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Next reverse the order of the input data x(2), x(3), x(4): 
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Careful examination reveals that the above matrix product conforms to the definition of dis- 


crete circular convolution: 

N-l 

y(k) = x(n) * h(n) = ^ x(n)h(k - n) (9) 

n c 0 

k = I, 2, * * < N 

Convolution is more easily understood from a graphical description as the mathematical process 
of equation (9) may not be clear. Figure 1 demonstrates graphically the convolution of x(n) and 
h(n), Both x(n) and h(n) are depicted as periodic discrete data scries of length four. For the case of 
k = 0,h(-n) is seen to be the mirror image of h(n), h(l - n) is simply h(-n) shifted right one sam- 
pling interval. y(l) may then be computed as the sum of products of x and h as shown. The other 
y(k) terms may be similarly found us demonstrated in the figure. 

Equation (8) describes the convolution: 

|x(i),X(2),X(4),X(3)| = | x(l), x(3), x(4), x(2) j * { w 1 , w 2 , w 4 , w 3 [ (10) 

The graphical details of this convolution may be studied in Figure 2. Note that upon convolving the 
two data sets of equation (10), the elements X(l), X(2), X(4), X(3) respectively are computed. Sig- 
nificantly, the elements X(k) of the DFT are now to be computed from a convolution operation. 
This transition was accomplished by a mapping of the matrix indices and is always possible for N 
equal to a prime or a prime power. For a mathematical description of the mapping process see 
Kolba and Parks (Ref. 9). 

FAST DISCRETE CIRCULAR CONVOLUTION 

Winograd (Ref. 6 and 7) has demonstrated an operational advantage to computing the DFT by 
changing to a discrete circular convolution operation. He presents an algorithm for performing 
short length cyclic convolution in a minimum number of multiplies. The concept employs polyno- 
mial multiplication modulo a third polynomial. 
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y(1) = h 0 x 0 + h 3 x 1 +h 2 x 2 + h 1 x. 


V(2) = h 1 x 0 + h 0 x 1 +h 3 x 2 + h 2 x : 


y(3) = h 2 x 0 + h 1 x 1 + h 0 x 2 + h 3 x 


y(4) = h 3 x 0 + h 2 x 1 + h 1 x 2 + h 0 x 


Figure 1. Graphical Convolution of the Periodic Data Sets x(n) and h(n) 



X(1 ) = w'xd) + w 3 x(3) + w 4 x(4) + w 2 x(2) 


X{2) = w 2 x(1) +w 1 x(3> + w 3 x{4) + w 4 x(2) 


X(4) = w 4 x(1 ) + w 2 x(3) + w 1 x(4) + w 3 x(2) 


X(3) = w 3 x{1) + w 4 x{3) + w 2 x(4) + w 1 x(2) 


Figure 2. Graphical Convolution of |X(1), X(2), X(4), X(3) } 
x(3),x(4), x(2) \ * jw 1 , w 2 , w 4 , w 3 ( 



The cyclic convolution h(n) * x(n) can be evaluated from the N coefficients of: 


Y(z) = U(z) • X(z) mod (z N - I) 


where 

N-t 

H(z) = h(.z k = h 0 + h j z* + h 2 z 2 + . . . + h N _ |Z N " 1 (1 1) 

k°0 
N-l 

X(z) ” ^ ^ X k Z k = X 0 + X j Z 1 + X 2 Z 2 + X N _ | z N " 1 

k°o 

Winograd states that the minimum number of required multiplies is equal to 2N - K where K is the 
number of irreducible polynomials into which z 14 - 1 may be factored, i.e., 

z N - 1 = n Qj(z) 02) 

M 

To demonstrate this theorem consider the convolution |h 0 , hj , h 2 , h 3 } * { xq> x ]» x 2 , x 3 | as 
presented earlier in Figure 1. 

y(l) = x 0 h 0 + xih 3 +x 2 h 2 + x 3 h, 
y(2) = x 0 h 1 +x | h 0 + x 2 h 3 + x 3 Ji 2 

1 1 j 

y(3)“ x 0 h 2 +X|hj + x 2 ir 0 4 *X 3 .h 3 
y(4) = x 0 h 3 +xjh 2 +x 2 hj +x 3 h 0 

The y(k) may be evaluated from the polynomial multiplication: 

Y|(z) = (x 0 + x l z + x 2 z 2 +x 3 z 3 ) • (h 0 + hjz + h 2 z 2 + h 3 z 3 ) 

= (x 0 h 0 ) + (x 0 li, + x,h 0 )z + (x 0 h 2 + x,h 1 + x 2 h 0 )z 2 + (x 0 h 3 + Xjh 2 +x 2 h, +x 3 h 0 )z 3 

+ (x,h 3 +x 2 h 2 +x 3 h,)z 4 +(x 2 h 3 +x 3 h 2 )z 5 + (x 3 h 3 )z 6 

- y 0 + yj z + y 2 z 2 + y 3 + y ^ + Yu 2 - 5 + Y6 z6 ( 14 ) 

It is now required to find Y(z) = Yj(z) mod (z 4 - 1). 
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The notion of modulo arithmetic is from the concept of congruency. Two integers a and b 
arc said to be congruent mod M if 

a “ b + kM (15) 

where k is an integer and M is the modulus. Tin’s may be written as: 

a = b mod (M) 

The integer b may be found as the remainder of the quotient a/M as may be easily demonstrated. 
Prom equation ( 1 5) we have 

a - kM ~ b 

and hence: 

k 

m nr 

kM 

b 

gives b as a remaindei. 

Intending this concept to polynomials we compute Y(z) = Y, (z) mod (z 4 - 1) as the 
remainder ui the synthetic division Yj(z)/(z 4 - 1). This results in 

Y(z) = (y 0 + y 4 ) + (y , + y 5 )z + (y 2 + y 6 )z. 2 + y 3 z 3 (16) 

Note that there arc N = 4 coefficients resulting from 

Y(z) = H(z) • X(z) mod (z 4 - 1) 

and that upon close observation they are found to be y(l), y(2), y(3), and y(4) of equation (13), i.e„ 

yo + y4 = x o h o +x i h 3 +x 2 h 2 + x 3 !i i = y(0 
yj +ys = x o ll i +x ] li 0 + x 2 h 3 + x 3 ii 2 = y(2) 
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y 2 +y6 = x 0 h 2 +xjh, +x 2 h 0 +X3I13 = y(3) 
y 3 = x 0 h 3 + Xjh 2 + x 2 hj + x 3 h 0 = y(4) 


The above example has demonstrated the logistics of performing the multiplication of two 
polynomials modulo another polynomial and specifically how the convolution results are recovered. 
This process may be accomplished in a more computationally efficient manner using the polynomial 
version of the Chinese remainder theorem (Ref. 9), We may first evaluate H(z) and X(z) modulo 
the irreducible polynomial factors of z N - 1, i.c., Qj(z) and then find Y(z) as 


where 


Y(z) = 


K 

X) Y I (Z)S 1 (Z) 
1=1 ' 


mod (z N - 1) 


(17) 


Yj(z) = H|(z)X|(z) mod (Q,(z)) 

Hj(z) - H(z) mod (Qj(z)) 

Xj(z) ® X(z) mod (Qj(z)) 

Sj(z)=l mod (Qj(z)) i = 1 1 2, • • • , K 

To demonstrate the operations described by equation (17) we may continue with the convolu- 
tion example |h 0) hj, h 2 . h 3 | + | x 0> Xj , x 2 , x 3 } as presented earlier. For this case N = 4, hence 

H(z) = h 0 + h j z + h 2 z 2 + h 3 z 3 
X(z) = x Q + x ( z + x 2 z 2 + x 3 z 3 

as before. 

Tlie irreducible real factors of z 4 - 1 are 

(z 4 - l) = (z + l)(z- l)(z 2 + 1) 

Note that the number of irreducible factors K = 3, and hence the number of required multiplies to 
compute the convolution is 2N - K = 5. 
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The Xj(z) and Hj(z) polynomials are determined next. Xj(z) is found by: 

X j (z) = Xg + x j z + x 2 z 2 + x 3 z 3 mod (z + I ) 

Computationally Xj(z) is computed as the remainder in X(z)/(z + 1). Similarly X 2 (z) and X 3 (z) 
may be found and \vc have: 

Xj (z) = X(z) mod (z + 1) = x 0 - Xj + x 2 ~ x 3 = xj 

X 2 (z) = X(z) mod (z - 1 ) = x Q + x j + x 2 + x 3 = Xq (18) 

X 3 (z) = X(z) mod (z 2 + 1 ) = (x j - x 3 )z + (x 0 - x 2 ) = + x pz 

Since H(z) lias the same form as X(z), the Mj(z) polynomials will have the same form as the 
Xj(z) polynomials: 

H, (z) « H(z) mod (z + 1) = h 0 - h j + h 2 - h 3 = hj 

H 2 (z) = H(z) mod (z - 1) = h 0 + lij + h 2 + li 3 = h^ (19) 

H 3 (z) = H(z) mod (z 2 + 1) = (h j - h 3 )z + (hg - h 2 ) = + h pz 

Wc may now find the Yj(z) polynomials as 

Y j (z) = H j (z)X j (z) mod (z + 1 ) « hj x 0 ! = yj 

Y 2 (z) = I-I 2 (z)X 2 (z) mod (z- 1) = hg Xq = (20) 

Y 3 (z) = H 3 (z)X 3 (z) mod (z 2 + 1) 

= hpxpz 2 + (ii^xp +hpx^)z + hgXg mod (z 2 + I) 

= (llQ X 3 + ll [ Xg )z + (llg 3 X 0 3 - IlpXp) 

-y 0 3 +y/z 
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The Yj(z) polynomials may be written as 


YjOO- ho x o = yJ 

Y 2 (z) = h 0 2 x 0 2 = y 0 2 (21) 

Yj(z) - (Iiq Xg - hpxp) + [(Hq - 1>1 3 )( X ? ‘ x o) + h O x O +ll i 3 x 1 3 ]z 
= y 0 3 +y, 3 z 

Next let 

m I =h 0 x 0 m 4 = ^ 0*0 

m 2 = hg Xg m 5 = hj 3 Xj 3 (22) 

m 3 = (Iiq - h p)(x p - x 0 3 ) 

These are the live multiplies that will be required to complete the convolution in this example. The 
Yj(z) polynomials may be rewritten in the form 

Y| (z) - nij 

Y 2 (z) * m 2 (23) 

Y 3 (z) - (m 4 - m 5 ) + (in 3 + ni 4 + m 5 )z 

We need yet to express Y(z) in terms of the Yj(z) factors as per equation (16): 

Y(z) = Y, (z)Sj(z) + Y 2 (z)S 2 (z) + Y 3 (z)S 3 (z) mod (z 4 - 1) 

To do so, however, we must first determine the S- (z) polynomials. These are found in accordance 
with equation (16) to be: 

S ] (z) = - Va(7? - Z 2 + Z - 1 ) 

S 2 (z) = !4(z 3 + z 2 + z + 1) (24) 

S 3 (z) = J / t (z 4 - 2z 2 + 1) 

11 



Those may be checked by verifying that S,(z) = l mo( | (q^)) !1s 

H4(z 3 -z 2 + z- !)= i mod (z+ 1) 

Va(7? + z 2 + 2 + 1)= i 111 oil (z - l) 

‘A(z 4 - 2z 2 + 1) = l mod (z 2 + )) 

Now Y(z) is found to be: 

Y(z) = m | [~A(z 3 - z 2 + z - l)J 
+ m 2 ['/i(z 3 +z 2 + z+ I)] 

+ [(m 4 - m 5 ) + (m 3 + m 4 + m 5 )z] [J4(z 4 - 2z 2 + 1)] mod (z 4 - I) 

4Y(z) « (m, + m 2 + ni 4 - m 5 ) 

+ (~ni | + m 2 + m 3 + m 4 + m 5 )z 
+ (m , + m 2 - 2m 4 + 2m 5 )z 2 
+ (-ni, + m 2 - 2m 3 - 2m 4 2m 5 )z 3 

+ (ni 4 - m 5 )z 4 

+ (m 3 +m 4 +m 5 )z 5 mod (z 4 - ]) 

4Y(z) = ( m , +m 2 +2m 4 - 2m 5 ) 

+ (-m, + ni 2 + 2m 3 + 2m 4 + 2m 5 )z 
+ (m , + ni 2 - 2ni 4 + 2m 5 )z 2 

+ (-m l + ni 2 ~ 2ni 3 " 2 »i 4 ~ 2m 5 )z 3 


(25) 


(26) 


(27) 


Recalling the method or determining one polynomial modulo another, equation (27) was found as 
the remainder in thesynthetic division 4Y(*)/ft< - I). As with equation (id), the coefllcients of a 
are the elements ol tile convolution example, i.e. : 
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( 28 ) 


yC 1) = !4(ni| + m 2 + 2m 4 - 2m s ) = x 0 h 0 + Xjh 3 + x 2 h 2 + x 3 hj 

y(2) = ‘/i(-mi + m 2 + 2m 3 + 2111 ^ + 2m s ) = x 0 h| + Xjhg +x 2 h 3 +x 3 h 2 
y(3) = ‘A ( nij + m 2 - 2m 4 + 2m 5 ) = x 0 h 2 + X|hj + x 2 h 0 + x 3 h 3 

y(4) = */»(— in j + m 2 - 2m 3 - 2m 4 - 2m s ) = x 0 h 3 + Xjh 2 + x 2 hj + x 3 h 0 

This Inst result may be verified by performing the indicated multiplications and additions. 

COMPUTING THE DISCRETE FOURIER TRANSFORM WITH FAST 
CONVOLUTION TECHNIQUES 

The fast convolution algorithm may now be applied to the task of computing the discrete 
Fourier transform. This is done by performing the indicated convolution in equation (8) to find 
X(l), X(2), X(4), and X(3). In this application h 0 , h, , h 2 , and h 3 are w 1 , w 2 , w 4 , and w 3 , respec- 
tively, as seen by comparing Figures 1 and 2. Similarly, x 0 , Xj , x 2 , and x 3 are x(l), x(3), x(4), and 
x(2), the input data vector. 

The terms hj , h 0 2 , h 0 3 , hp, xj, x Q 2 , x 0 3 and xp of equations (18) and (19) may be found as: 

hg = w 1 - w 2 + w 4 - w 3 = 2.236 
!lg = w 1 + w 2 + w 4 + w 3 = 1 .0 
hg 3 = w J - w 4 =-j 1.902 
h 3 “ w 2 - w 3 = -j 1.176 
Xg 1 - x(l) - x(3) + x(4) - x(2) 

Xq = x(l) + x(3) + x(4) + x(2) 

Xg = x(l) - x(4) 
x 3 = x(3) - x(2) 

From these results the nij terms in equation (22) may be computed and subsequently the con- 
volution results of equation (28) found. This results in the intermediate transform values X(k). The 
DFT is then calculated from equation (5). Tire process is now complete and the results for an N = 5 
computational algorithm are summarized in Table 1 [from Kolba and Parks (Ref. 9) and Winograd 
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Table 1 


Computational Algorithm for an N = 5 DFT Using Fast Discrete Convolution 
[from Kolbu and Parks (Ref, 9) and Winograd (Ref, 7)] 


a t = x(l) + x(4) 


a s = a 2 + a 4 

a 2 = x(l)- x(4) 


u 6 = a 1 - a 3 

a 3 = x(2) + x(3) 


a 7 = aj + a 3 

a 4 = x(2) - x(3) 


a 8 = x( 0 ) + ay 

m | = 0,951 a 5 


c, = x( 0 ) - m 5 

m 2 = 1.539 a 2 


c 2 = C[ +m 4 

m 3 = 0,363 a 4 


c 3 = c, -m 4 

m 4 = 0.559 a 6 


c 4 = m, - m 3 

m 5 = Va a ? 

X(0) = a 8 
X(l) = c 2 -jc 4 
X(2) — c 3 - jc 3 
X(3) “ c 3 + jc 3 
X(4) = c 2 +jc 4 

c 5 “ 1,1 2 “ m l 


(Ref. 7)]. The number of calculations required is observed to be 17 additions and 5 multiplications 
(the multiplication by Va may be accomplished by two word shifts on some FORTRAN compilers). 
The x(n) input data vector may be complex in which case these are complex adds and multiplies. 

Computational algorithms for other short length transforms may be found in Winograd (Ref. 7), 
Silverman (Ref. 8 and 10), and Kolba and Parks (Ref. 9). Tiie number of computations required for 
these several short-length WFT algorithms is compared with radix-2 FFT requirements in Table 2. 
The computational advantage of the WFT is clearly seen. 
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Tabic 2 


Number of Calculations for Short-Length WFT and FFT 


N 

WFT 

FFT 

Multiplies 

Adds 

MultipMes 

Adds 

2 

0 

2 

1 

2 

3 

2* 

6 



4 

0 

8 

4 

8 

5 

5* 

17 



7 

8 

36 



8 

2 

26 

24 

48 

9 

10* 

49 



16 

10 

68 

32 

64 


♦The number of multiplies may be reduced by using word shifts, 


For large N the WFT algorithm for fast convolution gets out of control. Furthermore, the 
number of additions becomes excessive and the WFT loses its computational advantage. Conse- 
quently it is necessary to use a combination of short-length transforms to realize the speed advantage 
of the WFT for transform lengths of practical importance. 

LONG-LENGTH TRANSFORMS 

Winograd Fourier transforms of practical lengths in the several thousands may be computed as 
a combination of short-length transforms. This is accomplished by converting an N - Mj M 2 * • * Mg 
(where the Mj arc mutually prime integers) length transform into £ shorter transforms of lengths M; 
for i = 1 , 2, • • • £. This is equivalent to a mapping from one to £ dimensions. 

The DFT of equation (1) 

N-l 

x (k) = x(n)W nk 

n=0 
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incorporates an index n which orders the input data vector and an index k ordering the transform 
output results, A mapping from one to two dimensions requires that each index map to two 
indices, i.e,, 


n^(n,,n 2 ) 

k - * (kj , k 2 ) 


With this mapping the DFT may be written as 


I ' k 2^ = 23 ( 23 X(n I ' n 2 J W "h 2 W M > 1 


k, = 0, 1, • • • • M, - I (29) 
k-, = 0, J 


where 


W Mi = o 




n=m,m 2 


and Mj and M 2 arc relatively prime. 

The two-dimensional DFT of equation (29) is found by firsv computing Mj transforms of 
length M 2 

Mi-l 

y('ij , k 2 ) = x(n, , I1 2 ) wj, 2 ^ 2 
H2=0 

and then Mo transforms of length Mj 

M j -1 

X( k i , k 2 )= 23 y(n l’k 2 )w M J kl 

ni=0 
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The above mapping follows the method of Kolbtt and Parks (Ref. 9) and is referred to as a 
prime factor FFT algorithm. Winograd (Ref. 6 and 7) presents another method of structuring short- 
length transforms to accomplish the same end. This technique is referred to as the nested algorithm. 

From an implementation point of view the major consideration between the two algorithms is 
the amount of computational effort. The computational requirements of these algorithms for 
several values of N arc compared in Table 3. In general, the nested algorithm requires fewer multi- 
plies, while tiie prime factor algorithm requires substantially fewer adds. The total number of calcu- 
lations is substantially fewer for the prime factored algorithm. 


Table 3 

Comparison of Prime Factor and Nested Algorithms 


N 

Factors 

Prime Factor Algorithm 

Nested Algorithm 

Multiplies 

Adds 

Total 

Multiplies 

Adds 

Total 

252 

9 . 7.4 

1024 

6344 

7368 

848 

7128 

7976 

504 

9*7-8 

2300 

13948 

16248 

1704 

15516 

17220 

1260 

9*5*7*4 

7136 

40288 

47424 

5168 

50184 

55352 

2520 

9*5*7*8 

15532 

86876 

102408 

10344 

106667 

117011 


The particular computer available and the relative timings for floating-point multiplication and 
addition would dictate which algorithm would be more time efficient in any application. The com- 
puter used in this study was the IBM 360/91 at Goddard Space Flight Center. The 360/91 is a very 
fast system utilizing a high-speed “cashe-pipeline”. Floating-point multiplication and addition 
require essentially the same execution time. Hence, the nested algorithm was judged faster and 
selected for application on this system. 

APPLICATION TO APERTURE ANTENNAS 

As discussed earlier, the far-fieid radiation pattern of an aperture antenna can be determined to 
a good approximation as the DFT of the aperture field distribution. In this application it appears 
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desirable to replace the FFT, which is currently used, with tiic faster WFT algorithm. Perhaps the 
savings in epu run time would permit extension of the analysis program to two-dimensional geo- 
metries with several thousand data elements along each axis. At this time a reasonable estimate of the 
maximum input data vector size is about five thousand, This would allow analysis of a one meter 
antenna at 180 GHz with half wavelength sampling. 

A FORTRAN IV program to calculate the WFT was adapted to run on the IBM 360/91 com- 
puter and extended to handle data vectors up to 5040 points with double precision arithmetic (See 
Appendix), Transforms may be computed for any data length given by the product of four mutually 
prime numbers selected from the integers 2, 3, 4, 5, 7, 8, 9, and 16 (e.g., N~9X8X7X5= 2520), 
As a basis for comparing the performance of this WFT program with the conventional FFT, a 
benchmark aperture distribution of special interest was selected. The distribution specified is sine 
cn a pedestal with a 20 db edge taper and is representative of focal-point fed parabolic antennas 
exhibiting axial symmetry. This distribution is described by 0.0909 + 0.9091 sin ( j ^), m = 0, 

1 , 2, • • • * N - 1 ; where N is the number of samples in the aperture. The antenna diameter is speci- 
fied as 20 feet and the wavelength as 0.44973 feet. 

A 2520 and 5040 point WFT analysis were performed on this antenna and the results com- 
pared with a 4096 point FFT computation. Theoretically the WFT and FFT algorithms give 
exactly the same results so this comparison is intended as a verfication of the correctness of the 
WFT program. Figure 3 presents the 5040 point WFT and 4096 point FFT results. The two 
patterns are essentially identical except for minor differences in the higher order lobes. This differ- 
ence is on the order of half a db and is entirely the result of differences in precision in the two 
programs. The WFT was run in double precision (8 byte data length) and the FFT in quartic preci- 
sion (16 byte data length). 

A 2520 and 5040 point WFT arc compared in Figure 4. The results are very nearly identical. 

The only significant difference is the amount of detail produced. Understandably, the 5040 point 
transform yields twice as much detail as the 2520 point transform. 
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SOLID LIME-5040 POINT WFT 
SYMBOL LINE-4096 POINT FFT 

WAVELENGTH --= 0.45' : APERTURE = 20' : 20db EDGE TAPER 



Figure 3. Far-field Radiation Pattern for Benchmark Aperture by 5040 Point WFT and 4096 Point FFT 



SOLID LINE-5040 POINT WFT 
SYMBOL LINE-2520 POINT WFT 

WAVELENGTH = 0.45' : APERTURE = 20' : 20db EDGE TAPER 
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4. Far-fleld Radiation Pattern for Benchmark Aperture fay 5040 Point and 2520 Point WFT 



Figure 5 demonstrates the new algorithm’s ability to handle complex input data. A 360° linear 
phase gradient was imposed on the benchmark aperture, resulting In the expected squint. The differ- 
ence in results as computed by the FFT and WFT is again due to the greatly extended precision 
specified in the FFT program. 

The interval timer available on the 360/91 system library at GSFC was used to time the WFT 
and FFT subroutines. For this comparison the FFT subroutine was rewritten in double precision 
math so that both transform algorithms would be judged on the same basis. The interval timer is 
represented as accurate to the nearest 0.01 second, The results of this test arc presented in Table 4. 
WFT and FFT timings arc competed for values of N that are as numerically close as possible given 
the different constraints on N for the two algorithms. The adjusted time ratio is determined by 
scaling the time ratio by the ratio of the number of data samples for the two algorithms. 


Table 4 

WFT and FFT Timing Results 


Number of Data 
Samples 

CPU Time 
(Sec) 

Time 

Ratio 

Adjusted 
Time Ratio 

WFT 

FFT 

1008 


0.06 






7.2 

7.0 


1024 

0.43 




2048 

0.90 






3.8 

4.6 

2520 


0.24 




4096 

1.90 






4.0 

4.9 

5040 


0.48 




Clearly the WFT shows a minimum 4.6 to 1 speed advantage over the FFT for the range of N 
considered. This is considered a significant improvement in speed. On the basis of these timings it 
is estimated that two-dimensional WFT computations would require about one minute, ten minutes, 
and 40 minutes for the 1008, 2520, and 5040 point transforms, respectively. Using the conventional 

it 
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SOLID LIME— 5040 POINT WFT 
SYMBOL LINE-4096 POINT FFT 

BENCHMARK APERTURE/— 360 DEG. PHASE GRADIENT 



Figure 5. Far-field Radiation Pattern for Benchmark Aperture with Negative 360° Linear Phase 
Gradient by 5040 Point WFT and 4096 Point FFT 




FFT program, approximately seven minutes, 3 1 minutes, and 130 minutes would be required for 
the 1024, 2048, and 4096 point transforms, 

To achieve the large N capability of the WFT program it was necessary to nest the summations 
described by equation (29) three and four deep. The limits on these summations arc the prime 
composites which make up N (e.g., 14 = 2X5X7X9 = 630), The ordering of these factors greatly 
affects the run time for a given N. Minimum epu lime is realized when the composite transforms 
arc ordered in some optimum way, Table 5 presents the optimum ordering for N = 504, 1008, 
2520, and 5040. Any departure from this ordering will result in an increase in run time by as much 
as 100 percent. 


Table 5 


Optimum Ordering of Composite Transforms 
for Large N WFT 


N 

Order 

504 

8*9*7* 1 

1008 

7* 16*9*1 

2520 

7*9*5*8 

5040 

16*5*7*9 


OTHER ALGORITHMS 

Other possibilities exist for fast DFT computations, Morris (Ref. 1 1) reports that a radix-4 
FFT algorithm outperforms the WFT on some computer systems (DEC PDP-1 1/55 and the IBM- 
370/ 168). His work indicates that in this comparison the WFT execution times were about 20 to 
40 percent longer. 

Reed and Truong (Ref. 12) have suggested that replacing the convolution operation in equa- 
tion (8) with a complex integer transform operation will result in fewer multiplications than either 
the FFT or the WFT. Several papers by these authors have demonstrated the feasibility of perform 
ing cyclic convolution by a combination of two N X N integer transforms, a multiplication of two 
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1 X N matrices, and an N X N inverse integer transform. The transforms and the inverse transform 
are performed with word shifts and integer additions. This results m a rather significant reduction 
in the total number of multiplications and converts them from floating-point to integer operations. 
However, the number of required additions is approximately tripled; consequently, the total num- 
ber of mathematical operations is increased. On the basis of these observations it appears that the 
use of integer transforms offers no advantage over either the FFT or the WFT for systems like tiie 
IBM 360/91 where floating-point multiplication and addition require about the same time. A clear 
advantage can be demonstrated for micro-processor based systems for which the ratio of multipli- 
cation to addition times may be several orders of magnitude. 

CONCLUSION 

Use of the WFT algorithm in antenna analysis appears to be a very successful application. The 
radix-2 FFT as used in computing far-ficld radiation patterns of aperture antennas may be replaced 
by the WFT with no degradation in performance and with a considerable improvement in speed. 
Over the range of N from about 1000 to 5000 the WFT demonstrated a minimum 4.6 to 1 speed 
advantage over the presently used FFT. This is sufficient to make two-dimensional transforms up 
to N = 2520 feasible within reasonable computer run time limitations. 

As new algorithms and transforms are introduced into the study of antennas, more powerful 
analysis and design techniques will become available to the design engineer. 
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APPENDIX 


FORTRAN PROGRAM FOR COMPUTING THE FAR-FIELD RADIATION PATTERN 

USING THE WFT ALGORITHM 


This section presents two FORTRAN subroutines. The first (GOODFT) uses the WFT algo- 
rithm to compute the DFT for complex input data up to length N = 5040. The second (PATOUT) 
computes the far-field radiation pattern from the transform results, The WFT subroutine is an en- 
hancement of a program contributed by Dean Kolba (Ref. 9) from Rice University. An N = 16 

WFT algorithm was added to the original program structure to extend the maximum range from 
* * 

‘ H ... 

N = 2520 to N = 5040. In addition, the program was rewritten in double precision, 

The length of the DFT, N, must be a product of no more than four mutually prime factors 

chosen from the integers 2, 3, 4, 5, 7, 8, 9, and 16, These factors are named Ml , M2, M3, and M4. 

If not ail four factors are used the unused factors are set equal to 1, The factors of one must be last 

in the sequence of M’s in the program. The other I/O variables used in the subroutine are; 


where 


NFT = number of nonunity factors 
KOUT - output indexing constant 

= K1 +K2 + K3 + K4 (mod N) 


K1 = M2*M3*M4 

or 

= 0 for Ml = 1 

K2*M1*M3*M4 

or 

= 0 for M2 = I 

K3 = Ml *M2*M4 

or 

= 0 for M3 = 1 

IC4 = Ml *M2‘M3 

or 

= 0 for M4 = 1 


XR(N) = r A al part of input data 
XI(N) = imaginary part of input data 
A(N) = real part of transform results 
B(N) = imaginary part of transform results 
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To illustrate the above, consider the case N = 5 *3 *2* 1 =30. The above input variables are: 


N = 30 


NFT = 3 


MI = 5 

K1 =6 

M2 = 3 

o 

a 

CJ 

' M3 = 2 

K3=T5 

M4 = ’l 

K4 = 0 


KO(JT = 3 l 'Cmod 30) = I ’ ' 

The ordering of the M’s is not significant with respect to the quality of the transform results, 
but does have a very important affect on the run-time for the subroutine. Table 5 presented the 
optimum ordering for the four cases N = 504, 1 008, 2520, and 5040. Departure from these orders 
may increase the cpu timings by more than 100 percent. 

The DFT is cyclic in nature and hence to compute the transform of a non-cyclic, finite data 
set jx(n)| it is necessary to append a relatively large number of zeros to both ends of the aperture 
distribution, About 75 to 90 percent of the input data should be made up of these embedded 
zeros. This will establish an adequate ground plane about the aperture and reduce the effects of 
folding or aliasing. 

The subroutine PATOUT requires as input the length of the DFT (N), the transform results 
from GOODFT (A&B), and the sampling interval (T) in wavelengths. To prevent aliasing T < A/2. 
The output of the subroutine is the magnitude (in db) and the phase (in degrees) of the antenna 
gain as a function of polar angle starting at a line drawn broadside to the antenna and through its 
center. For this statement to be true it is necessary that the input aperture distribution to 
GOODFT be specified according to the same geometry. The aperture data must be input starting 
at the antenna midpoint and proceeding to the edge. The zeros are imbedded next, followed by the 
other half of the aperture data starting at the edge and ending at the center. 
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SUBROUTINE GOODFT ( XR , X I , N, Ml , M2 » M3 , M4 , NFT , KOUT ♦ A , 6 ) 

C THE SUBROUTINE GOODFT COMPUTES A LENGTH N OFT OF THE INPUT DATA WHICH IS IN 
C TWO VECTORS, XR THE REAL PART AND XI THE IMAGINARY PART. BOTH XR AND XI ARE 
C LENGTH N VECTORS. THE LENGTH OF THE DFT * N, MUST BE A PRODUCT OF AT MOST 
C FOUR MUTUALLY PRIME FACTORS. THE POSSIBLE FACTORS ARE 2 , 3 , 4, 5 , 7 . 8 ,9 AND 16. 

C THESE FACTORS ARE Ml, M2, M3, AND M4. IF THE FOUR FACTORS ARE NOT ALL USED, 
C THE UNUSED FACTORS ARE SET EWUAL TO 1. FOR EXAMPLE WITH N=30, WE HAVE 
C MU5, M2 = 3, M3 = 2 » AND M4=l. THE FACTORS OF ONE MUST BE THE LAST OF THE M*$. 
C THF NUMBER OF NONUNITY FACTORS IS NFT. KOUT IS AN OUTPUT INDEXING CONSTANT 
C WHICH IS PRECOMPUTED. KOUT= ( K 1+R2+K3+K4 )MOO N WHERE K 1=M2#M3*M4 » 

C K2=M1*M3*M4, K3=Ml*M2#M4, K4=M1*M2*M3, AND K2=0 IF M2»l* K3=0 IF M3=l, AND 
C K4=0 IF M4= 1 . FOR EXAMPLE ,• N"30, Kl=6, K2 = 10, K3=15, K4=0 AND KOUT = 31MOO 30 
C =1. THF TRANSFORMED RESULTS ARE STORED IN TWO LENGTH N VECTORS, A AND R. A 
C CONTAINS THE REAL PART ANn R CONTAINS THE IMAGINARY PART OF THE -RESULTS. 
IMPLICIT RE AL*B (A-H.O-Z) 

DIMENSION XR I 5040 ) ,X 1(50*0) ,UR ( 14 ) ,U I ( 16 ) , I ( 16 I , A( 5040) ,B( 5040) 

RFAL*8 MR1 , MR2 , MR3 , MR4 , MR5 , MR6, MR7, MR8, MR9* MR 10 , MR1 1 , MR 1 2, MR 1 3 
REALMS MR14,MR15,MR16.MR17,MR18,MR19.MR20,MR21,MR22,MR23,MR24 
REALMS MR25,MR26,MR27,MR28,MR29,MR30,MI1,MI2,MI3,MI4,KI5,MI6,MI7 
REALMS MI 8, MI 9, MI 10, MI 1 1 ♦ MI 12 ,M 1 13 ♦ M 1 14 ,M 1 15, MI 16, MI 17, MI 18, MI 19 
REAL*8 MI20,MI21,MI22.MI23,MI24,MI25,MI26,MI27,MI28,MI29,MI30 

nf=nft 

c oroer factors for transforms of length Ml 

MM 1 *M 1 
MM2=M2 
MM3=M3 
MM4=M4 
GO TO 20 

10 GO T0( 12, 13, 14) ,NF 

C ORDER FACTORS FOR TRANSFORMS OF LENGTH M2 

12 MM1=M2 
MM2=M1 
MM3=M3 
MM4=M4 
GO TO 20 

C ORDER FACTORS FOR TRANSFORMS OF LENGTH M3 

13 MM1=M3 
MM2=M1 
MM3=M2 
MM4=M4 
GO TO 20 

C ORDER FACTORS FOR TRANSFORMS OF LENGTH M4 

14 MM1=M4 
MM2=M1 
MM3=M2 
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MK4=M3 

c indexing initialization for the transforms 

20 N2=0 
N3=0 
N4 = 0 

K 1 =MM?’t'MM3 , i , MM4 
K2=MM1*MM3*MM4 
K3=MM11‘MM2#MM4 
K4=MM1‘XMM2X«MM3 
I (1 )=0 

C INPUT INDEXING ALONG ONE DIMENSION 

21 00 22 0=2, MMl 

I ( J)=I< J-lI+Kl 

IF( 1 (U).LT.N) GO TO 22 

IIJ)=I(J)-N 

22 CONTINUE 

C TRANSFERRING DATA TO TEMPORARY VECTORS UR AND UI 

30 DO 31 0=1, MM1 
IJ=I ( J)+l 

UR ( 0 ) =XR ( 10) 

31 UI ( J ) -X I ( I J ) 

C TRANSFORM UR,UI 

GO TO ( 50,200,300,400,500,50,700, 800,900,50,50,50,50,50,50, 1600) , MM 
21 

C PL AC F RESULTS OF TRANSFORM BACK IN XR AND XI 

40 DO 41 J = 1 , MM 1 
I J = I ( J)+l 

XR ( I J ) =UR ( J ) 

41 Xi 1101=01(0) 

C TESTING FOR COMPLETION OF THIS FACTOR'S TRANSFORMS 
IF(W2,NE.MM2-1) GO TO 51 
N2=0 

I F( N3.NE .MM3-1 ) GO TO 52 
N3=0 

IF(M4.NE.MM4-1 ) GO TO 53 

50 NF=NF-1 

IF(NF.EO.O) GO TO 1000 
GO TO 10 

C INPUT INDEXING ALONG OTHER DIMENSIONS 

51 N2=N2+1 

DO 54 0=1, MM1 
I (0 )=I (0 >+K2 
IF( I(J).LT.N) GO TO 54 
I ( 0 ) = I ( 0 ) -N 
54 CONTINUE 
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nPT( '. T ’ \L I? AG A ^ 
01U' 1 * ArMll'Y 

Qp l 5 " f :- ALU 


r,n Tn 30 

62 N3=N3+1 

1(1 >=K3#N3+K4*N4 
IF( I ( n.LT.N) r,n TO 21 
1 ( 11 = 1 ( 1 ) — N 
Gn in 21 

53 NA=N4+1 
I I 1 )=K4*N4 

r,n rn ?i 

C 1 IN 5CW AM FI. ING TRANSFORM RESULTS 
1000 11=1 
J = 1 

GO TO 1001 

1002 IF(J.GT.N) GO TO 1003 
I I = I I +KDUT 

100A IF! I I .LE.N) GO TO 1001 
II=II-N 
Gn TO 1004 
1001 A ( J ) =XR ( I I ) 

13 ( J ) =-X 1(11) 

J=J + 1 

GO TO 1002 

C 2 POINT TRANSFORM 
200 URX=UR( 1 ) +UR ( 2 ) 

UIX=UI I1I+UK2I 
HR ( 2 ) =IJR ( 1 ) -UR ( 2 ) 

III (2 )=UI ( 1 )-UI( 2) 

UR ( 1 >=URX 
UI ( 1 > =U VX 
Gn Tn 40 

C 3 POINT TRANSFORM 
300 AR =UR ( 2 ) +UR ( 3 ) 

A I =U 1(2) +U 1(3) 

MR1=-1 .5D0*AR 
MI1=-1.5D0*AI 

MR2=O.R6602 5403RDO>X(UR( 2 >-UR(3) I 
MI2 = 0.nf>A025403800*<UI (2J-UI (3)1 
UR ( 1 ) =AR+UR ( 1 ) 

UI ( 1 )=A I +U I ( 1 ) 

MR 1=UR ( 1 1+MR1 
MI 1=UH 1 )+MI 1 
UR( 2) =MR1-MI? 

UI (2 >=MI1+MR2 
UR ( 3 ) =MR 1+M 1 2 
UI ( 3 ) =M I 1-MR2 
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on m 40 

C 4 (’DINT TRANSFORM 
AOO AR 1 =UR ( n+l)R(3) 

A 1 1 = U I (1)+UI (3) 

AR?=UR I 1)-UR(3) 

A I 2 = U l ( 1 > -01 I 3) 

AR3=UR ( 2 ) +UR ( 4 ) 

A I 3 = U I 12 j+UI I A) 
AR4=UR(2)-UR(4) 

AI 4=U 1 ( ? ) -U I ( 4 } 

UR t 1>=AR1+AR3 
UI{l)=An + A13 
UR ( 2 ) = AR2-A 14 
II I (?) “A I 2 +AR4 
UR ( 3 ) =AR1- AR3 
Uf ( 3 ) = A 1 1 -A 1 3 
UR(4)=AR2+AJ4 

UI ( 4 ) =A I 2-AR4 
GO TO 40 

C 5 POINT TRANSFORM 
500 AR 1=UR ( 2 ) +UR 1 5 ) 
AI1=UII?)+UI<5) 

AR2=UR I 2 )-URI 5 ) 

AI2=UI{ 2 I-UII 5) 

AR3=UR ( 3 ] +UR ( 4 ) 

AI3=UI (3)+UI(4) 
AR4=UR(3)-UR(4) 

AI4=UI (3 ) -U I (4) 

AR5=AR1+AR3 

A15=AI1+AI3 

MR 1=0. P5 10565 1630041 AR2+AR4) 
MI 1=0. 95 1 0565 163004 ( A I 2+A 14 ) 
MR 2=1 .538841 769004 AR 2 
M I 2=1 .63RR4 176900 4 Al? 

MR3 =0.3 632 71 2640004 AR4 

MI3=0.3632712640004AI4 

MR 4=0 , 5590169944004( AR1-AP3) 

MI 4=0. 5 5901699440041 A 1 1- A 1 3 ) 

MR5=-1. 25004 AR5 

MI 5=- 1 . 25D04A I 5 

UR ( 1)=UR( H + AR5 

1JII1 )=UI I 1 ) +A 15 

MR5 =UR ( 1I+MR5 

MI5=UI( 1 >+MI5 

AR1=MR5+MR4 



c 


A l 1=MI 5+M 14 
AR2=MR5-MR4 
A I 2=M I 5-M 1 4 
AR3=MRl-MR3 
A I 3=M I 1-M I 3 
AR4=MR1~MR2 
AI4=MI 1-MI2 
UR(?)=AR1~AI3 
>11 ( ? ) =A I 1 + AR3 
OR ( 3 ) =AR2+ A 14 
UI< 3 I =A 1 2-AR4 
IJR ( 4 1 = AR2- A I 4 


OElGTSAL 

. , , > • rt 

rv ± 


PAGE IS 

OVALITY 


UI ( 4 1 =A 1 2+AR4 
IJR ( 5 ) = AR 1 + A I 3 
U! ( 5 > =A I 1-AR3 
fiO TO 40 

7 POINT TRANSFORM 
700 ARl=UR<2)+UR(7J 
A 1 1 =U I (2 )+UI { 7) 

AR2=UR ( 2 1 -UR I 7 > 

A 1 2=U II 2 ) -U H 7 ) 

AR3 a UR ( 3 ) +UR 1 6 ) 

A 1 3=U I < 3 ) +U I ( A | 

AR4 = UR ( 3 I -UR ( A ) 

AI4=UI (3 ) -U I ( A) 
AR5=UR(4)+UR(5) 

A I 5-U I ( 4 ) +U I ( 5 ) 

AR A=UR ( 4 ) -UR ( 5 ) 

A I A=U I ( 4 1 -U I ( 5 ) 

AR7=AR1+AR3+AR5 

A I 7 = A 1 1 +A I 3+A I 5 

MRl=-l. IAAAAAAA7D0KAR7 

MI1=-1. 1A6AA6A6700*AI 7 

MR 2 =0.790156468800*1 AR1-AR5 > 


MI 2=0 . 7 90 1564 48800* ( A II** A 15) 
MR3=0. 05505426800* I AR5-AR3) 

MI 3=0. 055854260004) A 15-AI 3) 

MR 4=0. 734302 20 1D0*{AR3-ARI ) 

MI 4 = 0. 734 30220100 4 ( A I 3- A I 1) 

MR 5=0. 44095855200*1 AR2+AR4-ARA) 
MI 5 = 0. 44095855200*1 A I 2 + A I 4-A I A ) 
MRft=0, 34 08729310 0*( AR2+ARA) 

MI 6=0.3408729 3100*1 A I 2+ A I A) 


MR7 = -0.5339A936 100* ( - AR A-AR4 ) 
MI 7=-0. 53396936 IDO* ( -A 1 A- A 14) 
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MRrtnO. R7A8A229inO#l ARA-AR ^ ) 
Mlfl a n.87A8A2?91D0>!<( AIA-A12) 
DRIlhIJRI l)+AR7 
UI ( 1 ) =tJ I ( 1 ) +A 1 7 
AR 1 =UR I 1 1 +HR 1 

All =U I ( n+MI 1 

AR?=AR 1+MR2+MR3 
A J 2®A 1 1 +MJ 2+M 1 3 
AR3 =AR 1 -MR2-MRA 
A I 3 = A 1 1-M I 2-M J A 
AR A=AR 1-MR3+MRA 
A I A“A 1 1-MI 3+M I A 
AR5=MR8+MR6+MR7 
A I 5-M I 5+M I ft+M 1 7 
AR6-MR5-MR A-MRH 
A I I 5-M 1 fi-M I H 
AR7=MR9-MR7+MR8 
Af7»MlS-MI7+W|H 
UR ( ? ) aAR2-A I 5 
U! ( ? ) nA I 2+AR5 
IJR(3)=AR3-AI6 
IJI ( 3 ) =A I 3+AR6 
UR { A ) = ARA+ A 1 7 
UT (A ) =A r A-AR7 
(JR ( 8) =ARA-A 17 
UI ( 3 )=A I4+AR7 
IJH I ft ) = AR3 + A 1 6 
U! (M-AI3-AR6 
UR ( 7) =AR2+AI5 
UI I 7 ) =A 1 2-AR5 
on m ao 

C fi PDINT TRANSFORM 
Ron AR 1-UR ( ? ) -UR ( 8 ) 

A 1 1 =U H ? ) -U 1 1 8 ) 
AR2nUR(?)+UR(8) 

A I7-U I ( ? ) +U I (81 
AR3=UR (A1-URI6) 

A I 3-U I ( A ) -U I (61 
AR A=UR ( A ) +UR ( 6 I 
A I A = U I (AJ+UMft) 

AR 3 =UR ( 1)-UR(5) 

A I 5=U I ( 1 )-UI (5) 

AR 6=UR ( 1 ) +UR ( 5 I 
A 1 6=U I ( 1 )+UI I 5) 

AR7=UR 1 3 )-UR< 7t 



ORIGIN^ WjSBfc 


A 1 7«U t < 3 I -UI C 7 1 
ARfl=UR(3)-*UR(7) 

A! B=UI ( 3 ) +UI I 7 ) 

MRl=O.7O710A7Rl?0Q>M AR1+AR31 

MIl=0.7071067Rl?n0>MAl 1+ a I 3 > 

MR?cO. 70T10A7R'l2n0’«( AR2-ARA) 

Ml 2=0. 70710A7R1 200>> ( A 1 2-A 1 4 1 

MR3=AR?+AR4 

M I 3 = A 1 2+A I A 

MRA=ARMARfl 

M 1 4c A l A+ A I 8 

MR5=ARA-ARR 

MI 8 = A J A-A I 8 

MRA = AR — A R 3 

MI A=A I l-A 13 

MR7=AR5+MR2 

MI 7 = A I 5+M 1 2 

MRR=AR5-MR2 

MIR=A I 5-MJ2 

MR9=AR7+MR1 
MJ9=AI7+MIl 
MR I0=AR7“MP 1 
MI 10=AI7-MI1 
UR I 1)=MR4+MR3 

UI m=MIA+MI3 

UR ( ? ) =MR7-M I 9 
UK? )=MI7+MR9 
UR ( 3 ) =MR5-M 1 4 
UII 3 ) =M I 5+MRA 
UR ( A ) =MRR+M 1 10 
III < A ) =M I fl-MRIO 
UR( 3)=MRA-MR3 
UI ( 8 ) =M I A-M 1 3 
UR ( h ) =MR8-M 1 10 
UI ( A ) =M I H+MR10 
UR I 7 ! =MR*>+M I b 
U J ( 7 )=M | 5-MRA 
UR ( R ) =MR7 + M 1 9 
UI ( R ) =M t 7-MR9 
60 TO AO 

C 9 POINT TRANSFORM 
900 AR 1 =l)R ( 2 ) +UR ( 9 ) 

A 1 1=UI ( 7. )*t-UI I 9 ) 

AR?=UR( 2)-UR (9) 

At 2=U I (?) -U II 9 ) 



AR3=lJA(3|+URm 

M3=UI m+UI I R) 

AR4a|jR ( 3 ) -UR ( R ) 

AJA=UM 3MJIIA) 

AR5=UR(5)+UR(61 
A I 5=U I ( 5 I +IJ ! I M 
AR6=UR ! 5 )— UR ( 6 ) 

A I A=U I ( 5 I -LI I 1 A ) 

AR7oURtA)+UR(7) 

A 17MH (41+UI (7) 

ARR=UR ( A ) -UR ( 7 ) 

A I R = U I ( A ) -II I ( 7 ) 

AR=AR1+AR3+AR5 

A I =A I 1 +A I 3+A [ 5 

MRi=-n.5no*AR7 

Mll=-n.5no»Ai7 

MR?=n.RAAn?5AO3RD0>!>ARH 

MI2=0.HA60?5AO3Rno*AIR 

MR3=n. 1974654?n0*(-ARl+AR5) 

MI 3 = 17. 1974654?(7n*l-Al 1+AI51 
MRA=O.3Aft5790?n0«l AR1-AR3) 

MJ A = 0,SAR5790?fl0>K{ A I I - A I 3 > 

MR R =0, A 711 136170* I -AR3 + AR5 ) 

MJ 6=0. 371 111AI7n*I -A I3 + A1 5) 

MR6=0. 54 253 179170*1 AR2-AR6) 
MI6=O.64253l79D0*< A J2-AJ6) 

MR7=0. 1002567900*1 AR2+AR4 ) 

MI 7=0. 100255/900*1 AI2+AI4) 

MR fl=0. 44 22759700*1 -AR4-AR A', 

M 1 A=0 . 4422759700* I - A 1 4-A • 6 ) 

MR9=-1.500*AR 

M19=-l. 5I70*AI 

MRlOsO. A66O25403K0O*( AK2-AR4+AR6 ) 

MI 10 = 17. fl 66025403900*1 A I ?-A 14 + A I 6) 

AR1=UR( l)+MRl 

A II = U I ( n+Ml 1 

DR | 1 )=AR+AR7+UR< l I 

mu )=Ai+Ai7+uim 

AR=UR I 1 1+MR9 

A I =U I ( 11+MI9 

AR2=MR4— MRS 

A 1 ? = M 1 4-M I 5 

AR3=MR3+MR4 

A I 3®M 1 3+M f 4 

AR4=MR7-MRB 
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AI4=MI7-MIfl 

AR5=MRfc-MR7 

A t 5 = M I 6-M 1 7 

AR6=AR?~MR5-MR3+AR1 

AI6aAI2-MI5-Mt3+AIl 

AR7 S AR3+MR3+MR5+AR 1 

AI7=AI3+MI3+MI5+AI1 

ARR=-AR3-AR2+AR1 

Alft=-AI3“AI2+Al 1 

MR 1 =MRA-MRR 

Ml 1 = MI A-M 18 

MR3 =AR4+MR l+MRp 

MI3-AI4+MIUMI2 

MR4=AR5+MR 1-MR2 

Ml 4=A 1 9+M 1 1-M I 2 

MR5«AR5-AR4+MR2 

MI 5=A I 5-A I 4+M I 2 

UR I ? ) ~ AR6-M 1 3 

UI ( 2 ) -A 1 A+MR3 

UR ( 3 ) = AR7-M 14 

UI ( 3 ) =A I 7+MR4 

l)R { A ) = AR-M I 10 

UI ( 4 } s A I +MR 10 

UR { 3 ) = AR8-M 1 5 

tJ! (5 ) — A I 8+MR5 

UR ( 6> = AR8+M I 5 

UI ( A ) =A I A-MR5 

UR(7)-AR+MII0 

IJII 7 ) =A 1 -MR 10 

UR ( A > = AR7+M 14 

UII S ) =A I 7-MR4 

l)R ( 9 ) = AR6+M I 3 

1J1 < 9 ) =A I A-MR3 

GO TD 40 

C 1 A POINT TRANSFORM 
1 AOO AR 1 =UR ( 1 » +UR ( 9 I 
A 1 1 =U I ( 1 )+UI ( 9) 
AR2-UR(5)+UR( 13) 
AI2=UI (51+UM 13) 
AR3=UR ! 3 )+UR (11) 

A 1 3=U I ( 3 ) +U I ( 11) 
AR4=UR ( 3 I -UR ( 11) 

A 1 4=U I { 3 ) -U I ( 11) 
AR5-UR(7)+UR( 15 ) 

A I 5=U I ( 7 ) +Ul ( 15 ) 
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AR 6 =UR( 7 )-IJR< 15 l 
A I 6=U I ( 7 ) -U I ( 15 ) 
AR 7 =UR ( 2 ) +UR 110 ) 
A 17 =UI (2 l+UK 10 ) 
ARB=l)R I 2 ) -UR ( 10 ) 
A!fl=UJ 1 2 )-U I ( 10 ) 
AR 9 =UR I 4 ) fUR ( 1 ?,) 

A 1 9 = 1) I (4)+UI ( 12) 
AR10=UR(4) -UR( 12) 
A J 10 = U I (4 |-UI( 12) 
AR11=UR< 6)+UR< 14) 
A 11 1=1) I ( 6J+UI ( 14) 
AR 1 ?=UR ( 6 ) -UR ( 14) 
A I 12=UI (61-UI ( 14) 
AR 13=UR ( R ) +UR ( lft ) 
A J 13=11 1 ( 8 ) +U I (16) 
AR14=UR(R)-UR( 16) 
AI 14=U1 ( R )-UI ( 16) 
AR 1 5=AR 1 + AR2 
A 1 1 5 = A 1 1 +AI 2 
AR 1 6=AR3+AR5 
A 1 16 = A I 3 + A I 5 
AR 1 7=AR 15+ AR 16 
A 1 17 = A 1 15+A 1 16 
AR 1 R=AR7+AR 1 1 
A 1 1H = A I 7+AI 1 1 
AR 19=AR7-AR 1 1 
AI 19 = A! 7- A I 11 
A R 2 0*AR9+A R 1 3 
A I 20=A I 9+AI 13 
AR? 1=AR9-AR13 
A I 21 = A 1 9- A 1 1 3 
AR22=AR1R+AR?0 
A I 22 = A 1 1 R+A I ?0 
AR 23 = *RR+AR 14 
AI23=AIR+AI 14 
AR24=ARR-AR14 
A 1 24 = A I R-A I 14 
AR25=AR10+AR12 
A I ?5 = A I 1 0+A 1 1 ? 
AR?6=AR12-AR10 
A I 26 = A I 12-A 1 10 
AR 3 1 =(JR ( 1 ) -UR ( 9 ) 

A I 31 =U 1 tl)-UI ( 9 ) 
UR ( 1 ) =AR 1 7+ AR2 2 



OttlGltfAk 
x? pnou 


page is 

quality. 


Hill ) =A ! 1 74-A I 22 
HR < 9>=AR17-AR22 
III (9)=AI 1 7~A I 22 
AR29 = AR 1 5- AR 16 
A I 29 = A I 15 -A I LA 
AR30=AR 1-AR2 
A I 30 = A I 1-AI? 

MRl=0.7071067!U?nO*( AR19-AR21 ) 
MI 1=0.7071067*1200*1 A I 19-AI2I ) 
MR2=0.7071067R12DO*{ AR4-AR6 ) 
MI2=0.7n71067fil2n0*( AI4-AI6) 
MR3=0.3fi26H3A32400*( AR24+AR26) 
MI3=0.3A26a3432400*{ AI24+AI26) 
MR4=1 .3DA5A29A5nO*AR2A 
M 14=1. 3065 62 965 D0*AI?4 
MR5=-n.54119A100inO*AR?A 
M I5=-0. 5411961 001 DO *AI26 
AR3? = AR 1R-AR20 
AI32=-AI 1R+AI20 
AR33=AR3-AR5 
A I 33 = -A I 3 + A 1 5 
AR34=l)R ( 5 ) -HR ( 13) 

A 1 34=-U I ( 5)+IJl ( 13) 

MR6=n. 7071067M1200*! AR 19+AR21 ) 

M I 6=-0 . 7Q71D67H1 2D0*( A I 19+AI21 ) 

MR 7 =0 . 7071067R 12D0* ( AR4+ AR6 ) 

MI7=~0.7O7in67R12D0*( AI4+AI6) 

MRB=0.9?3R79532500*( AR23+AR25) 

MJfl=-0.9238795325D0*( AI23+AI25) 

MR 9 =-0.541 196100100*AR23 

M 19=0. 54 U96 100 100* A I 23 

MR 10=1 ,306562965D0*AR25 

MI 10=-1.306562965DO*AI25 

MR1 1=AR30+MR 1 

M 1 1 1 = A I 30+M 1 1 

MR 1 2 = AR30-MR 1 

MU2=AI30-MI1 

MR13=AR33+MR6 

M 1 1 3 = A I 3 3 +K 1 6 

MR 14=MR6-AR33 

MI 14=M1 6-AI33 

MR 15 = AR31 + MR? 

MI15=AI31+M12 
MR16=AR31-MR2 
M 1 1 6= A 1 3 1 -M 1 2 
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MR17=MR4-MR3 
M 1 1 7 = M 1 4 -M 1 3 
MR l ft =MR 5 -MR 3 
MI1R=MI5-MI3 
MR 1 9=MR 15 + MR 17 
MJ19sMI15+MI17 
MR20=MR 15-MR17 
MI2n=MJ 15-MI17 
MR 2 1 =MR 1 A+MR 1 R 
M 1 2 1 =M 1 1 6+M I 1 R 
MR 22=MR 1 A~MR 1R 
MI 2?=M I 1A-M1 1R 
MR23=AR34+MR7 
MI 23 =4 I 34+MJ7 
MR24 = 4 R34-MR7 
MI 24 = A I 34-M I 7 
MR25=MRR+MR9 
MI 24 = M I R+MI9 
MR2A=MRR-MR10 
M I ? A=M I fl~MI 10 
MR27=MR23+MR25 
Ml?7 = MI ? 3 +M I 2 5 
MR?R=MR23-MR25 
M I ?fl =M l 23-M 1 25 
MR29=MR24+MR26 
M I 29=M I 24+M I2A 
M R 3 0 »M R ? 4 — M R 2 4 
M{ 30 = MI 24-MI 25 
IJR ( 2 ) =MR 19+M 127 
IJI ( 2 ) =M 1 1 9+MR27 
UR (3 ) =MR 1 1+M 1 13 
UI ( 3 ) =M I 1 1+MR 1 3 
UR ( 4 ) =MR22-M 1 30 
UI (4 ) =M I 22-MR30 
UR ( 5 ) sAR29+A 132 
UI I 5 ) = A I 29+AR32 
UR ( A) =MRRi+MI29 
UI ( A ) =M I 21+MR29 
UR ( 7 ) =MR1 2+M 1 14 
UI (7 ) =MI 12+MR14 
lJR(ft)=MR20-M!2fl 
UI ( R ) =M I 20-MR 2R 
UR ( 10) =MR20+MI2R 
UI ( 10 ) =MI 20+MR2R 
IJR ( 11 ) =MR 1 2-M I 14 
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UK LI )=MI 1P-MR14 
UR ( 1?)sMR?!1-M1?<} 

IK ( 12 ) I Zl-MR?*} 

UR< 13 ) =AH?9-AI 3? 

UI (1 3 ) - A I 29-AR3? 

UR( 14)=MR2?+MF30 
If! ( 1 A 1 =M I 22+MR30 
UR ( 15) =MR1 l-MI 13 
UI ( L 5 ) =M 1 1 1 — M R X 3 
UR ( 1 A ) =MR19-M 1 27 
UI ( 1M sM F 19-MR27 
on in 40 
1003 RFTURN 
FWO 

SIJHR OUT I ME PATOUT (N, T I 
IMPLICIT REAI.AR (A~H,n~?) 

COMMON/ OFT /AI 5040) , R( 5040) 

WR1TF! 4,45 ) 

45 FORMAT! 1H1,43X, 25HFAK F IFLO ’ANTENNA PATTERN/ 7X , 5HANGLF , 4X , 7HMAG ( PR 
2) , 5X . 5HPHASE ) 

R1?»A{1 )*#2+B{ 1 )**? 

PFG=57 , ?95 7795 1 3 1 DO 
on 50 1=1, N 
S-l I~1 )/<N*T> 

IFIS.GT. l-.O) GO TO 40 
C? = A( I }***2+5I I )$*? 

CDR=io.ono«nLnGin(c?/ci?) 

PHA=D AT AN ( B ( I ) / A { I ) )*OEG 
IFI5m.LT. 0 . 0 ) PHA=PHA+ 150.0 DO 
ANG=PARSIN! S)*DFG 
WRITE! 4, 70 ) ANG, COH , PH A 
50 CONTINUE 
40 CONTINUE 

70 FORMAT ( 3 ( 3 X , FI 0 • 4 , F 1 0 . 4 , F 10 .4 ) ) 

RFTURN 

FND 


it 


41 



